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I did off-lattice 3D Monte Carlo simulations for polymer translocation through a 

c 

^ \ narrow pore at low external field, trying to be as close to a direct approach as possible. 

' The process was found non equilibrium globally, but dynamics of the monomers close 

to the pore (the fold) was found close to quasi equilibrium. I observed a tension (or 
pressure) buildup near the pore, and subdiffusion for the displacement of the middle 
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monomer in the pore, as reported in the literature. I did not find the distribution of 
the reaction coordinate after some time of free diffusion predicted by the fractional 



. diffusion equation. 
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I. INTRODUCTION 



Translocation of polymer molecules through nanopores is an important process in 
biology. 5 Examples are the transport of DNA and RNA molecules across pores of the nuclear 
membranes, injection of the virus genetic information, gene swapping etc. Translocation of 
single-stranded DNA/RNA through a pore forced by external voltage was demonstrated 
in vitro as a single molecule signal by Kasianowicz et alA. In such experiment voltage is 
applied across a membrane with a single pore in it. Due to it negative charge, the DNA 
or RNA molecule is forced to drift across (see the topical review by Meller~). While inside 
the pore the polymer blocks the ion current, thus increasing the resistance, which effect 
can be detected with a single event resolution. Such experiments showed sensitivity to the 
nucleotide composition, etc. 4 . It is now believed that by using the blockage current only, no 
sequencing can be achieved. Yet electrophoretic translocation through a hole combined with 
another method is still one of the prime candidates for fast DNA sequencing in the future 
(see the review by Zwolak 34 ). Maybe most experiments use a pore made by self assembly 
of a-hemolysin in lipid bilayer. Such pores are so narrow ( 1.4 nm diameter in the most 
narrow place), that only a single stranded DNA or RNA can pass. Typically the RNA/DNA 
is of the order of hundreds base units, the applied voltage is few hundreds mV (or less) , and 
the blockage times are hundreds of microseconds. More recently experiments with silicone 
membranes have been reported.- A big development in silicone nanopore production has 
been made in the last years. A pore of 1 nm diameter is reported, which similarly to the a- 
hemolysin can translocate only single stranded DNA-. Obviously wider pores can be created 
and used, but this is outside the emphasis of this paper. 

The translocation of a polymer through a narrow pore driven by electric field is also 
extensively simulated. Some authors use MD simulations with simplified (e. g. cylin- 
der like) poTe^-^^^^M^.^ some consider more detailed pore structure (a-hemolysin)i2>22>2!. 
Another approach is dynamic Monte Carlo simulations (only with simplified cylinder /square 
like p OTe ^£L2&2ZiM£l.£l£l. The last approach allows the case without external field to be 
simulated^^ 1 - 1 ^. No friction forces can be considered in a natural way in Monte Carlo 
simulations, hence the hydrodynamic interactions are neglected. Using some tricks that 
shorten the computer time, MD simulations on the low external field problem have been 
recently reported^ 1 ^ 
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The low external field behavior observed in simulations invoked some controversy. Con- 
sidering equilibrium entropic barrier one gets for the translocation time, reflecting boundary 
conditions at the beginning r = CN 2 , where N is the number of bases of the polymer (or 
other parameter proportional to the length of the chain). For ideal chain C = where D 
is ID diffusion constant. In the above mentioned result no friction is considered^* 8 - 1 ^. How 
noticed 8 , the non forced translocation time has a simple lower limit - the time the polymer 
diffuses its radius of gyration (or translocation with no barrier). That gives r > where 
D is the diffusion constant of the whole molecule. In Rouse model D ~ giving r > N 2u+1 . 
v is the Flory exponent, which in 3D is v = 0.5 for ideal chain, and v = 0.588 for self 
avoiding random walk, corresponding to polymer in a good solvent. The self avoiding case 
contradicts r ~ N 2 . 

Also, for long enough chains, the polymer shouldn't have time to equilibrate at each step 
of it translocation. Equilibration time goes as n 1+2u with the number of monomers outside on 
one of the sides of the membrane n. For ideal chain this gives N 2 , making equilibrium at each 
step possible if not probable, and N 2A8 in the self avoiding case, resulting in equilibration 
slower than translocation. 

Extrapolating the result of 2D Monte Carlo simulations, a r ~ N 1+2u big N behavior, 
good solvent was proposed^ (reflecting beginning, self avoiding chain, no external field). 
Later this result was supported in 3D by MD^^^. Based on 3D lattice^l^ 1 ^, or off 
lattice^ Monte Carlo much higher scaling exponents have been proposed for the above 
system, reaching up to^ r ~ N 2+u = jV 2 - 588 . Diffusion time scaling longer than squared 
length means subdiffusive behavior. 

I performed off-lattice Monte Carlo simulations equivalent to direct trial by transparent 
assumptions. The whole dynamics of the process was obtained this way as transition rates, 
and then r was calculated. I obtained for the above described system a smaller scaling 
exponent limited only by r > N 2u+1 in Rouse model. 

II. THE MODEL 

I regard the mean translocation (or first passage) time r of polymer through a hole in a flat 
membrane. I start with the first monomer passed, it cannot exit back (reflecting boundary 
conditions), and finish when one monomer is left not translocated. The membrane is thin 
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only ~ 1 monomer is in the pore, the others are outside. The hole is cylindrical and narrow, 
so the chain cannot enter folded, and in the theoretical fit, the monomer inside is assumed 
without degrees of freedom. 

The polymer consists on N beads. I use dynamic off-lattice Monte Carlo model close to 
the bead-spring model^^ 1 ^ 1 ^, which gives Rouse dynamics^. Instead using springs, I just 
give some allowance for the distance between the nearest neighboring monomers. For the 
self avoidance I assume the bead is a hard sphere with radius that prevents self crossing. 
I pick a monomer at random and attempt to move it ±0.5 at each x, y, and z direction 
simultaneously with uniform distribution of the probability. The move is rejected if it places 
the beads to close to (far from) one another or to the membrane, or accepted otherwise. 
This is equivalent to or oo potential, or there are no attractive or repulsive interactions. 
N such attempts are the unit of time (one MC step). No external field is applied - all the 
motion is completely stochastic by nature. 

I regard self avoiding and phantom chains. For the self avoiding case the distance between 
two nearest neighbors is between rmin = 0.5, and rmax = 0.7. With respect to not nearest 
neighbors, and the membrane the bead is a hard sphere with radius 0.3, which makes the 
self crossing of the chain practically impossible. The membrane is 0.11 in thickness, so no 
additional care is needed for the polymer not to drift with no bead in the pore. A thicker 
membrane slows down the translocation. The pore is cylindrical with radius 0.4 -just a little 
bit bigger than the bead. 

For the phantom polymers only nearest neighbors distance is controlled, self crossing is 
allowed. The distance between two nearest neighbors is again between rmin = 0.5, and 
rmax = 0.7. The membrane is 0.71 in thickness for the same reasons as above. The pore is 
cylindrical with radius 0.1. 

I checked the two models by using a polymer of n monomers with the first one fixed on 
the membrane surface, r connects the first and the last bead. First I compared the mean of 
the component of r in direction perpendicular to the membrane. Similarly to the bulk case 
this should give scaling as n u . I also checked the equilibration time of such chain defined 
as the correlation coefficient ffTOl) between f and it initial position as a function of the time. 
As mentioned one expects scaling as n l+2v . I did this for chains of n = 294 and n = 6 
monomers for each model. For the phantom case I got v ~ 0.5 within our accuracy from the 
both tests, as expected. For the self avoiding case I got v ~ 0.6 from the distance test, and 
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n 2 ' 3 as an approximation of n l+2v from the correlation test. Thus I overestimated a little 
bit the theoretical expectation of v = 0.588. Note that for similar polymer lengths similar 
scaling is reported for dynamical exponent^ 2 - as a reasonably close to the theoretical one 
(which should be observed when n — > oo). The obtained exponents proof that the polymer 
is effectively in a good solvent. 

The number of translocated monomers is the reaction coordinate k. I actually simulate 
the effective rates for the polymer to jump step monomer forwards or backwards. Then one 
can calculate the translocation time r from ([1]) and (|2J). I use step = 6 monomers (unless 
otherwise specified). This is much (about a order of magnitude for the regarded systems) 
faster than to wait for a direct escape starting from the beginning. I start from k = k and 
simulate till the chain reaches both k = n + step and k = k — step at least once. Than I 
start from k = h + step. The program keeps the configuration of the polymer when first 
reached translocation to k = k + step and use it as the initial state when starting from the 
new k = k + step state. I use (jSJ) to calculate the effective rates u e J l fecttve and w^ ectwe . 

I did 500 or more trials for every simulated point shown, unless otherwise mentioned. 

I define k = k ± step or k when the corresponding monomer is at the middle along it 
translocation in the pore. In the off-lattice MC monomers jump randomly, which introduce 
some error in defining the position. Reducing the relative error, a bigger step makes the 
result more accurate, but demands longer computer time. 

Please note, that in this work I regard the mean translocation time. This is the case for 
most theoretical works and simulations. In experimental works the most probable passage 
time is measured. For theoretical works dealing with the shape of the translocation time 
distribution (not zero external field) see^^ 3 -. 



III. THEORETICAL APPROACH 



I use the formula for the mean time r to first reach b or e. Let's b is before e, and we 
start from s between b and e. For discrete ID walk— 

HM jJ ,e +1 )^, l -M, e + 1) |l|li nj) 

\k=b k=b i=k+l j=k+l J / 

ts — l s—2 -. s— 1 i \ 

E^E^Eny (i) 
k=b k=b i=k+l j=k+l J / 
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where 

7T(0 — 1, s, e + 1) = r (2) 

1 Z^k=bllj=b Uj 



is the probability starting at s to reach e + 1 before ever reaching 6—1. Above is the 
rate per unit time to jump forwards from k to k + 1. it; & is the same but backwards from 
to k - 1. 

We also need the ansatz (see Fig. 1) 

effective _ n (^m, S m , e m ) 

m ~ r(b s e ) 

w effective _ \ ~ ^fei s m; e m) 

&m Sm+li b m Sm—1 (3) 

the last is consistent with the formula for mean first passage time, or if one uses u^^ ecUve , 
and w^ ecti ^ instead of u k , and w k in Eqs.flU), one gets the same time. Note that the set of 
u^ ectme , and w e J l ^ ecttve should be less in number than the one of Uk, and Wk in order (jHJ) to 
be needed. 

To estimate the rates theoretically I use: 



D (F{k) -F{k + 1)\ D 



A 2 V 2&gT ; A 2 



D f F(k)-F(k-l) \ D Z k -! 
Wk = A^ 6XP ( 2k^f =^ X '— (4) 



where D is the ID diffusion constant if the distance traveled for a jump between 2 neighbor- 
ing states is A. F(k) is the free energy at state k. The factor of 2 in denominator insures 
the detailed balance between the neighboring states. 

One may assume a smooth continuous along k potential and put additional states between 
the original ones, shortening A. One sees from ([3]) that D should be kept constant to 
reproduce u e J L ^ ectlve and yj^ff ectlve between the original states the same as the original rates 
if F is constant, and A — > 0. This way one obtains from (j3j) and ([I]) the continuous limit 
for the case we start from state 0, with reflecting boundary conditions there (wq = 0): 

N k 

r(0,0,"U- = /*«P g$) /*5T XP ("SB <5 » 
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For polymer of n segments attached at one of its end to a flat surface the partition 
function is 29 : 



Z n ~ n 



7/— 1 



Hence the entropic potential of our system has the form 
F(k) 



1 - 7/) ln(fc) + (1 - 7/) ln(iV - M - k) + const 



(6) 



(?) 



In the last equation M is the number of the monomers inside the pore, which is assumed 
without degrees of freedom, k is the number of translocated monomers or our reaction 
coordinate, and: 

( 1 

phantom chain 



< 1 

2" 



I 1 = \ 0.69 - self avoiding chain 



(8) 



1; 



rigid rod 



According to (j4j) and ([7]) our system have the symmetry: 



Uk — WN-M-k 
Wk = UN-M-k 



(9) 



How we know^&^ 
for ideal chain. 



with (El) gives r(0, 0,N — M) 



UIQ=0 



C{N - Mf with C = ^ 



IV. HOW EQUILIBRIUM IS THE TRANSLOCATION PROCESS? 

I first use the data for phantom chain. How noted, the translocation process for them 
is expected to be closer to quasi equilibrium one. All the simulated data shown are fitted 
to the theory by one constant. I calculate and simulate the first exit time for the polymer 
starting with one base passed, which cannot reenter back - w\ — (reflecting boundary 
conditions), and finish at one monomer left - r(0, 1, N — 1 — M) W1=Q . In all the data M = 1. 
Note that in this case the forward reaching probability tt — 1 in (JT]). 

In Fig. 2 is shown a fit of the effective forward rate u e J l ^ ectlve for ideal polymer to jump 
step = 6 monomers forwards, N = 303. The theoretical values are calculated using (j3J), 
([7]) and then ([3]) to get the effective value for step = 6 monomers. One sees almost perfect 
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agreement with quasi equilibrium theory. The simulated data for r(0, 1, N — 1 — M) Wl=0 as 
function of the number of the monomers in the chain - iV are shown in the Fig. 3 for the 
case of ideal chain. Since the rates agree very well with theory, I simulate only half of them, 
than use their symmetry (J9j) to find the rest (and r(0, 1, N — 1 — M) Wl=Q ). The slope goes 
down with N, reaching N 2,2 . In this graph one sees again a very good agreement with the 
quasi equilibrium discrete process theory, given by (jlj), (JTj), and (pQ). I believe for bigger N 
the exponent will drop to 2. 

The next question is if these data correspond to quasi equilibrium process. To estimate 
a lower limit of the equilibration time of a chain of n monomers anchored to the pore at one 
end, I use the correlation coefficient c(n,t): 



In the above definition r is the vector from the end of the polymer fixed in the membrane 
(the pore) to it free end, and ( ) means averaging over ensemble or over time which I 
assume equivalent. c(n, t) is correlation coefficient of the direction of r , hence it should 
vanish the same as or faster than the real memory of the whole chain. I did simulations for 
anchored chains of n = 294 and n = 6 monomers for both self avoiding and phantom models. 
The observed behavior is c(n, t) = exp(-u^j) where for n —>■ oo one should have u = l + 2u, 
and u is specific for each model. As noted for the phantom case I really have the big n 
value L) — 2, while for the self avoiding model u = 2.3. Similar approach and results with 
slightly different correlation coefficient, and using lattice dynamics are reported^ for the 
same system. For the phantom case I show the times t e (n) and tio(n), (where c(n,t e ) = -, 
and c(n, t w ) = in Fig. 4. Together are the mean time for the polymer to move step = 6, 
monomers forwards or backwards which are directly simulated, as well the mean lifetime 
the polymer to move one monomer left or right (which come as from the theoretical 
values from (jlj), and ©, also used to get the theoretical fit in Fig. 2). One sees that (with 
probable exception of near the center), there was not been a real equilibrium. 

To see why in such case the data look like equilibrium again the phantom chain is used. 
In Fig. 5 effective forward rate u e J l ^ ectlve is given for the polymer to jump step = 6 monomers 
forwards, N = 303, and iV = 807. This is the very beginning of the process (k close to 
0). I also give the theoretical prediction for the rate with the same translocated number of 
polymers, but replacing the rest (non translocated end) of the chain with much shorter or 




(10) 
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longer one (total number becomes N = 40, or N = 4000). One sees, that this doesn't change 
much the translocation rate. Or, if the monomers on one of the sides of the membrane have 
a counter part on the other of few tens of monomers or of few hundreds, this doesn't change 
much the dynamics around the pore. It worths to mention, that (HI) with ([71) saturates to 
limiting values for the beginning of the process (or by ([9]) the end of the process) as N grows: 

~ 1—7/ 
D ( k \~ 

\k+ l) 

~ 1 — 7/ 

D ( k \ 2 
A2 \k-l) 

small, N ^ 00 (11) 

So, I believe that the process is equilibrium for chains of few tens of monomers. In case of 
bigger chains, the ones close to the pore do not feel the rest of the chain more than few tens 
of monomers apart. In other words the rates are close to ones if everything greater than few 
tens monomers each side of the pore is cut. How suggested^ 1 ^ there is a partial equilibrium 
for the few tens of monomers close to the pore. 

In Fig. 5 one also sees that u^ ect%ve is consistently lower than the theoretical prediction. 
I interpret this fact as that there was a local equilibrium, and the rest of the monomers (not 
being equilibrated), just pull slightly backwards. 



u k = 

w k = 
k - 



V. DOES THE SELF AVOIDANCE KEEPS THE THEORY? 

Next I present data for self avoiding chain. Again they all are fitted with one constant. 
Due to the different lengths of the pores (thicknesses of the membranes), this constant is 
incompatible to the one for phantom chain. 

In Fig. 6 are given the the effective forward rate u f ^[^ ectlve for the polymer to jump 
step = 6 monomers forwards for N = 75, and the same for step = 30 monomers forwards 
for TV = 123. One may see that for the larger chain, there is some deviation from the 
theory. That is why I do not use ([9]), but rather simulate all the rates needed to calculate 
the translocation time. For the same reason, I also give in Fig. 7 times obtained by direct 
trials. One may observe the nearly perfect match between the theory, and the two ways of 
simulation. The translocation time r(0, 1, iV — 1 — M) Wi= q, M = 1 ends up scaling as N 2 . I 
conclude, that the hypothesis for local equilibrium next to the pore is still valid. 
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What violates slightly the theory is that the number of translocated monomers is not 
a perfect reaction coordinate. It doesn't give all the information for the state of system, 
namely that the polymer is stretched when it just advances (see the next section). 

VI. SHORT TIME DYNAMICS OF SELF AVOIDING POLYMER 

The procedure described above (but with step = 2 which shortens the computation time) 
is used to thread a self avoiding polymer exactly to the position where the middle monomer 
is in the pore. Then I set the time to zero, and observe the further dynamics. The entropic 
potential is then zero, and doesn't affect the process. 

It is believed that for small enough times the dynamics is subdiffusive. The data for 
SD 2 of the displacement of the middle monomer along the pore really scale as t 0J for the 
range shown in Fig. 8, N = 127. Interestingly the SD 2 for the reaction coordinate k scales 
linearly for this range. The displacement of the middle monomer along the pore also scales 
differently than the change of k. The author explains this with the stretching of the chain. 
Unfortunately very short time intervals are hardly accessible. In the off lattice Monte Carlo 
monomers jump to random positions, and for small times the error in defining k will become 
compatible to it value. 

Within the accuracy, there is no deviation from Gaussian distribution of the reaction 
coordinate k at some fixed time after starting from the middle (for the range shown in 
Fig. 8). The last is in agreement with the lattice simulations described^ where they have 
observed minute deviation far from the mean value for small times. This is also in agreement 
with the simulations-- , for phantom polymer in ID. Deviation from Gaussian distribution in 
the middle - a non smooth cusp is predicted by the fractional diffusion equation, describing 
subdiffusion with diverging mean waiting times between the consecutive moves. So, my 
results are in parallel with the statement^ that for the polymer translocation we may have 
a different type of subdiffusion. 

The average shift of the reaction coordinate k is shown in Fig. 8. One sees that the 
average position consistently goes back. That is why we didn't get perfect agreement with 
the theory for N = 123. The reason is the accumulated tension in the chain. My results for 
the mean displacement of k due to the relaxation of the over stretching confirm the tension 
relaxation exponent a = introduced in^. The evolution of the SD(k) is completely 
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diffusive for the times shown in Fig. 8 within our accuracy. In parallel with what I said the 
dynamics is again local, SD of k is practically independent on iV - see Table 1. The data 
for the average shift k shown in Fig. 8 are less independent on N. 

The author speculates that such dependence of the dynamics from the history is the 
explanation of the subdiffusive behavior expected for longer chains. 

VII. WHAT IS THE LONG CHAIN TRANSLOCATION TIME SCALING? 

The last point is to check why the exponent is lower than the theoretical limit r > N 2u+1 
for the Rouse dynamics. I performed a direct simulations (200 trials) for the self avoiding 
polymer using pore much bigger than the polymer length (or practically absent). I did trials 
for N = 125 and N = 175. In agreement with our equilibration times (using the correlation 
coefficient (II 01) ) the time scales as iV 2 ' 3 . Comparing them with the self avoiding data for 
normal pore shown on Fig. 7 they are ~ 50 times smaller. I estimated, that following such 
behavior this time will surpass the one for the regular data for iV ~ 10 12 . In other words the 
data are too far from the point where the the scaling N l+2u will become important. Let's 
note that it is hard to extrapolate the same physics for such big N. On the other hand, all 
the simulation data (including the ones presented here ) start with larger exponent, which 
goes down with increasing the N. Here an exponent lower than the limit for the Rouse 
model is reached, which suggests that r ~ is the real exponent. Also, Rouse model (and 
hence MC) overestimates the hydraulic friction in D. For the real systems, D ~ -^7 (Zimm 
exponent) which doesn't contradict r ~ iV 2 . For arguments that the Zimm exponent may 
not be valid for our problem, see^ 3 ^. 
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time 


N = 127 


N = 179 


N = 299 


MC steps 








10 4 


0.80 


0.79 


0.76 


10 5 


2.6 


2.6 


2.7 


10 6 


8.3 


8.2 


7.8 



Table 1. The SD(k) of a self avoiding polymer with time, due only to diffusion. SD(k) = 
at t = 0. 
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Figure Captions: 

Fig 1. Using the effective rates u e ^ ectlve and yj e ^ ecUve turns the system into equivalent 
one, consisting of smaller number coarse grained states. 

Fig 2. Simulated and theoretical rates u e ^ ectlve , step = 6 as function of the number of 
the monomers passed. The theoretical fit is done using quasi equilibrium approach. 

Fig 3. Translocation time r as a function of the length of the chain N. One monomer 
stays in the pore, process starts with one base passed, which cannot reenter back - w\ — 
(reflecting boundary conditions), and finish at one monomer left. The theoretical fit is done 
using quasi equilibrium theory for discrete states. 

Fig 4. Equilibration time at one of the sides of the pore (thin lines) and the time to pass 
step monomers during spontaneous translocation (thick lines). Upper thin lines is the time 
the correlation coefficient ( flOl) to vanish down to c(n, t) = the lower thin line - the same 
for c(n, t) — -. Upper thick line is the average time the polymer of N = 303 needs to move 
step = ±6 monomers which is simulated to get u e ^ ecttve and w e ^ ecUve , lower thick line - the 
same for step = 1. The diamonds are for TV = 807 and step = ±6. 

Fig 5. First few rates u e ^ ecUve , step = 6 for big phantom chains. The lines are calculated, 
as the original polymer is replaced with one having much shorter non translocated tail (total 
length N = 40, upper curve), or much longer non translocated tail (total length TV = 4000, 
lower curve). Quasi equilibrium discrete theory is used. 

Fig 6. Simulated and theoretical rates u e ^ ectlve as function of the number of the monomers 
passed. The theoretical fit is done using quasi equilibrium approach. In the case of N = 75, 
step = 6 and the fit is good. In the case N = 123 step = 30, and we have performed 
only 200 trials. In the last data one sees a consistent deviation from the theory, due to the 
accumulation of tension in the chain. That is the reason that step = 6 cannot be used. 

Fig 7. Translocation time r as a function of the length of the chain N. One monomer 
stays in the pore, process starts with one base passed, which cannot reenter back - w\ — 
(reflecting boundary conditions), and finish at one monomer is left. The theoretical fit is 



16 



done using quasi equilibrium theory for discrete states. 3 of the times are simulated using 
average of 200 direct trials. 

Fig 8. The polymer was threaded in using the procedure described in the text, step = 2, 
equivalent to spontaneous translocation. When the middle is reached, the further dynamics 
is observed, which is due entirely to the accumulated tension. The line corresponds to t a 
with tension relaxation exponent a = as predicted by-22. 
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Figure 2. Kotsev 
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Figure 7. Kotsev 
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Figure 8. Kotsev 



